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Abstract 

In this article we present a non-linear dynamic programming algorithm for the com- 
putation of forward rates within the maximum smoothness framework. The algorithm 
implements the forward rate positivity constraint for a one-parametric family of smooth- 
ness measures and it handles price spreads in the constraining dataset. We investigate the 
outcome of the algorithm using the Swedish Bond market showing examples where the 
absence of the positive constraint leads to negative interest rates. Furthermore we inves- 
tigate the predictive accuracy of the algorithm as we move along the family of smoothness 
measures. Among other things we observe that the inclusion of spreads not only improves 
the smoothness of forward curves but also significantly reduces the predictive error. 
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1 Introduction 



Within the financial industry forward rate curves play a central role in fixed-income derivative 
pricing and risk management. Notwithstanding, such curves are not empirical or directly 
measurable objects but rather useful abstract concepts from where observed prices can be 
derived. Moreover, given a finite set of market prices, we can construct in general an infinite 
number of compatible forward rate curves 1 . To avoid such ambiguity several approaches have 
been proposed in the literature trying to capture a "reasonable" or "natural" functional form 
within the set of compatible possibilities. 

In historical order the first kind of methods proposed to solve this problem make use 
of the so-called parametric approach. In this approach a particular functional form for the 
forward rate curve is assumed leaving a certain number of free parameters to be fixed from 
the calculation of a given set of quoted prices. An extensive literature exists advocating for 
this approach. We can cite as examples the works of McCulloc 1_, Vasicek and Fong 0, 
Chambers, Carleton and Waldman pj. Shea [I], Nelson and Siegel |5| and more recently the 
works of Svensson 6 , Fisher, Nychka and Zervos [7j and Waggoner (Hj. In most of these works 
we notice the privileged role played by polynomial and exponential splines as the preferred 
functional forms for the forward rate curves. 

The second kind of methods has been termed in the literature as non-parametric or 
maximum-smoothness approach. Here instead of advocating for an a priori functional form 
for the forward rate a given measure of smoothness is chosen and then the forward rate curve 
is obtained as the one maximizing this measure subject to the constraints imposed by market 
prices. Examples where these methods have been investigated include the works of Adams 
and Van Deventer j^j, Delbaen and Lorimier ^U], Kian Guan Lim, Qin Xiao et al [111 I12j . 
Frishling and Yamamura |13j and Yekutieli ^3]. In these works three different smoothing 
measures have been proposed. We also have the works of Forsgren |15j and Kwon |16j 
that generalize these methodologies and clarify the connection between splines and certain 
smoothness measures. Finally we point out the work of Wets, Bianchi and Yang ^2] that 
can be located somewhere in-between both approaches since here the number of functional 
parameters is finite (albeit arbitrarily large) and the functional behavior is restricted to a 
subfamily of C 2 curves. 

The purpose of this article is two-fold: Firstly we want to present an efficient maximum- 
smoothing algorithm that handles the presence of spreads and implements the positivity 
constraint. Secondly we want to investigate the predictive power of a linear combination 
of two quadratic measures, namely the one proposed by Delbaen et al. and Frishling et al. 

1 We assume here that prices do not allow for arbitrage. If this is not the case the compatibility condition 
has to be relaxed. 
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|10| 113] and the one by Adams and van Deventer [2] . Here it is worth remarking that once 
the compatibility with market prices is fulfilled the only guiding principle that should be 
taken as definition of "reasonable" or "natural" is the predictive power and not other ad-hoc 
criteria. 

In this article we will only use as constraining data coupon bearing bonds. The inclusion 
of treasury bills, zero coupon bonds or bill futures is straightforward and amounts to adding 
the corresponding linear constraints. Since our objective in this article is focusing on an algo- 
rithm dealing with non-linear constraints and inequality constraints (spreads and positivity 
constraint) we have not included such data. 

With these objectives in mind we organize the article as follows: 

In section |2] we present the objective function that we will use throughout the article 
and we establish the basic notation. In this section we present a sketch of the complete 
algorithm leaving the details for the appendices. In section we present the results of the 
article including examples where the absence of the positivity constraint or the adequate 
spreads leads to negative rates. Here we present also a study of the predictive power of the 
one-parametric family of smoothness measures that include as extreme cases the measures 
used by Delbaen et al., Frishling et al. and Adams and Van Deventer. Finally in section 0] 
we present the conclusions. 



2 The algorithm 

A bond, j, is an instrument that gives future coupons, Cij, at time stages ,i = 1, • • • , n,j— 1 

and a final payment, c nj j = Nj 2 . The bond price, Pj, can be determined from the discrete 

(i) 

forward rate curve, f r , r = 1, . . . , Rn J , as follows 



i=l 



\ 



r=l 



Mr , (1) 



where £ r is the length of the time period between time stage r and r+1 (in our implementation 

we have used £ r = 1 day). 

The objective function, or smoothing measure 3 , is defined as a linear combination of the 

one used by Delbaen et al. and Frishling et al. (DF) jlOl lloj and the one used by Adams and 

2 By final payment Nj we mean the complete last cash flow of bond j, typically that includes a principal 
plus a last coupon. 

3 Note that since we do not include a global minus sign we have to perform a minimization and not a 
maximization. With this sign, that is the one used in the literature, "rugosity" measure would be a more 
appropriate name. 
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Van Deventer (AD) 
n-l 



n— 1 / ^ \ 2 

Jr+l — J"? 



n-l 



(Cr-1 + £r 



r+1 



/r /r /r— ] 



£r-l 



£r, (2) 



r=l v / r= 2 

The first term in Eq.(J2J) (DF) is a discrete approximation of the integral of the square of the 
first derivative of / and the second term (AD) is a discrete approximation of the integral of 
the square of the second derivative of /. This objective function is to be minimized subject 
to the consistency constraints 



(3) 

St > 0, P) < Pi < P% (4) 

where pj = In (Pj/Nj) is to be determined along with / and where we have used the definitions 



r=l 
h 



■■ ^2 ay exp 



i=i 



^ ^ r ' P) := 111 ( a 



# := In (ag,) 



a:. 



0.7 



pa 

3 



a 



L, (5) 

with Pj, Pj - the respective bid and ask prices of bond j (J = 1, • • • , to) 4 . Note that the 
equality constraint given by Eq.Q is just Eq.Q rewritten taking logarithms and using the 
definitions (J5J). Eq.@ introduces two inequality constraints. The first one is the positivity 
constraint over the forward rate curve and the second one is the requirement that the single 
price given by Eq.lfT]) must lie in-between the bid and ask prices. 

We define the spread of bond j as the quantity p^ — p^. We take the largest time to 
maturity in Eq.pjl equal to the largest time to maturity in the constraining dataset, namely 



n := max 

j 



The constraints reflecting bond prices (pQ) have been rewritten in a way such that they 
become linear when no coupons are present (vj = 1 for a zero-coupon bond j). Constraints 
given by Eqs.(j3J) and @ are moved to the objective function defining 



Z:=W+J2 X i 

i=i 



Pj + f r & ~ ln ( v i 



r=l 



-p 1" (fr)-p £ (in [p 3 - pfj + ln (p<j - Pj 
r=l j=l 



(6) 



Note that this optimization problem is non-convex and therefore several local minima may exist. 
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with the Lagrange multipliers A,-, j = 1, ... ,m and the logarithmic barriers with parameters 
p > and p, > (in the solution procedure we take p — > 0, p — > 0). The use of log barriers 
to deal with inequality constraints is a standard methodology in interior point methods for 
optimization problems JH]- An explanation of this methodology adapted to our problem is 
given in appendix A. Briefly the minimization algorithm is structured as follows: 

step 0: Initialize log barriers with coefficients fj\°\ pP\ and set (/, p) = (f^°\p^) 

with the seed satisfying the inequality constraints @. Let k = 0. 

step 1: Make a second order approximation of Z at (/ (see appendix A), 

step 2: Determine the newton step (/ 

[fc+l] )j5 [fe+i]) using 

dynamic programming (see appendix B). 
step 3: Modify (/ [fe+1] , p [k+1] ) to get a solution (f [k+1] , p [k+1] ) that satisfies the inequality 
constraints Update log barriers (p^ — ► and — > as — > oo) and 
calculate the values of W and of constraints ©• Check if a termination criterion 
is satisfied, otherwise let k = k + 1 and go to step 1 (see appendix C). (7) 

Computing times involved in step 2 are summarized in subsection B.l. The solution typically 
stabilizes in approximately 6 iterations as can be seen in Fig.Q. On a Pentium 4, 2.4 Ghz 
computer the algorithm coded in C++ takes around 1/4 sec. to compute a forward rate 
curve like any of the ones seen in Fig.©. 

3 Consistency and predictability 

In this section we present some examples of the behavior of algorithm Q and we investigate 
the predictive power of the resulting forward rate curves. In Fig.(j2J) we present a series of 
forward rate curves calculated using DF and AD smoothing measures. There we can see 
that for both measures the resulting curves share some similar traits like the positions of 
most peaks and dips. Clear differences between both sets of curves are found at their end- 
points and in their behavior in the presence of high spreads. In the set of curves obtained 
from the DF measure we have vanishing first derivative at the end-points and curves that 
tend to constants for high spreads. For the AD measure we have vanishing second and third 
derivatives at end-points p2| and curves for high enough spreads given by straight lines. From 
the financial point of view these features are, in principle, just different aesthetic possibilities. 
In order to choose a particular measure the guiding principles should be, in the first place, 
the fulfillment of the consistency constraints given by Eqs.(jnHU) an d after this is guaranteed 
the predictive performance. 

Let us start now with the analysis of consistency. As can be found in ^] if we do not 
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Sweden, July 2001 




Figure 1: Plot of the objective function W (see Eq. (J2J) vs. Newton iterations (starting from 
the seeds: = 0.04/year, p^ = {p a - + p b j)/2) . The plots correspond to the forward rate 
curves of Fig.©. The converge behavior is affected by the termination criterion dictated by 
Eq. (|29|) and by the initial values of the log barriers dictated by Eq. (|30|) . Note the exponential 
convergence of the algorithm in the first steps. Note also that on Friday 06 the minimal 
solution for AD measure is, for all practical purposes, a straight line. 
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Figure 2: Forward rate curves obtained from the prices of 11 Swedish government bonds in a 
two- week period of July 2001. The curves in plot (a) are calculated using DF measure (7 = 1 
year 3 , (f = 0) while the ones in plot (b) correspond to the AD measure (7 = 0, 99 = 1 year 5 ). 
In both plots we have used a tolerance spread of 1% (pj — P b j = 0.01). 
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consider the positivity constraint the local minima of objective © are given by exponential 
splines with exponents iyT/V or polynomial splines of order 2 or 4 when ip or 7 are re- 
spectively zero. The main problem with these exponential or polynomial splines is that there 
is no warranty that they fulfill the positivity constraint. Negative rates are not admissible 
in the absence of arbitrage opportunities and the risk of obtaining this unwanted feature 
is illustrated in Fig.©. In this figure we have concentrated on the Swedish bond data on 
Monday, July 09, 2001. There we have tested three spread patterns for both DF and AD 
measures with and without the positivity constraint. From these plots it is evident that the 
inclusion of spreads in the calculation of the forward rates is a necessary ingredient that can 
have a major impact in the resulting functional behavior. 

Once we have an algorithm that insures the consistency of forward rates we can concen- 
trate on the predictive accuracy of different measures. However, before starting the analysis 
of this issue let us make a brief digression to comment a point regarding measure (|2"|). If we 
want to have both 7 and ip different from zero and we want to compare the effects of each 
term it is important to realize that DF and AD measures scale differently under changes of 
time units. In other words, 7 and ip have different units. A practical way to define their units 
is to consider the objective W as an adimensional quantity. By doing so and remembering 
that / has units of inverse time, it is immediate to obtain that 7 has units of time 3 and cp 
units of time 5 . The importance of keeping this in mind becomes apparent in results like the 
ones presented in Figs.© and ©. 

Figs.© and © summarize our results regarding the predictive accuracy of the algorithm 
as a function y^/c/?. There it is clear that the characteristic time span where DF and AD 
compete is not the day or the century, but clearly the year. To construct these figures 
we have calculated the forward rate curves for different values of y7/V when one bond is 
removed from the constraining dataset. The price of this missing bond is used afterwards 
as a benchmark to test the accuracy of the resulting curves. Since we are interested in 
the statistical performance we have done such comparison for 335 consecutive trading days 
starting on Wednesday, November 08, 2000 and ending on Thursday, March 07, 2002. 

We are also interested in studying the impact of spreads in the constraining dataset over 
the predictive accuracy. Therefore we present our results for three spread patterns, namely 
constant spreads of 0%, 0.5% and 1% in the constraining dataset. Fig.© concentrates on the 
predictive accuracy for the first 9 bonds of Table © and Fig.© presents the same analysis 
for the remaining 2 bonds of Table ©. These last 2 bonds are the ones with the larger 
maturities in the complete dataset. In particular for the last one with the largest maturity 
we have to decide upon the methodology to extrapolate the forward curve outside the range 
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of the constraining dataset 5 . Therefore in Fig. ([5]) we present the prediction accuracy using 
constant extrapolation from the last maturity in the constraining dataset and 'W-generated" 
extrapolation that consists in utilizing the VF-optimal forward rate curve even outside the 
range of the constraining dataset. 

4 Results and conclusions 

In this article we have presented a non-linear dynamic programming algorithm designed for 
the calculation of positive definite forward rate curves using data with or without spreads. 
We have included multiple details of the algorithm aiming at practitioners not familiar with 
the techniques of dynamic programming. We have illustrated the results of this algorithms 
using the Swedish bond data for a one-parametric family of smoothness measures. 
The results and conclusions are the following: 

• The proposed algorithm calculates forward rate curves in real time and admits, with- 
out time or complexity penalizations, the use of any non-linear local objective function. 
Since it also handles non-linear constraints it is possible to include within the constrain- 
ing dataset any derivative products with prices bearing some dependence on forward 
rates. 

• This is the first algorithm proposed in the literature that implements the positivity 
constraint in the maximum smoothness framework. To the knowledge of the authors 
the only other work that implements such constraint outside this framework is the one 
of Wets et al. JZj. The proposal in has the advantage of using simple linear 
programming but do not consider the presence of spreads minimizing instead the sum 
over the modulus of the difference between calculated prices and market prices. 

• For the objective functions and constraining datasets like the ones we have used or more 

generally for the ones studied in the algorithm proposed in that reference offers 

better computing times at the expense of ignoring the positivity constraint. Essentially 

the complexity in |16j is c 3 and in ours is n c 2 where n is the number of time steps and 

c the number of constraints. For that reason, when this class of objective functions and 

constraining datasets are used, a well coded algorithm might try first the proposal in 

|16j (improving its treatment of the spreads using e.g. log-barriers) and later, only if 

the result is not positive definite, use our approach. 

5 By the range of the constraining dataset we mean the range in time to maturity that goes from present 
to the last maturity within the dataset. 
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Figure 3: Forward rate curves obtained from the prices of 11 Swedish government bonds on 
Monday, July 09, 2001. The vertical lines indicate cash flows. Highest lines indicate the 11 
cash flows at maturities (normalized to 1) and the smaller ones coupon cash flows. In these 
plots we show the effect of the positivity constraint with three spread patterns both for the 
DF (in (al), (bl) and (cl)) and the AD (in (a2), (b2) and (c2)) measures. In plots (cl) and 
(c2) we see that a spread of 1% is enough to generate positive and sensible curves for both 
measures (the positivity constraint has no effect in this case). In plots (al) and (a2) we take 
null spreads observing that for both measures large negative rates are obtained (conspicuously 
for AD measure). Also in these plots we see that even though the positivity constraint can 
be fulfilled the resulting curves have large oscillations. Negative rates or large oscillations in 
the positively constrained curves can be related to price patterns that are close to violate 
absence of arbitrage (the precise nature of this relation will be studied elsewhere). Plots (bl) 
and (b2) explores this fact keeping null spreads in two particular bonds (bonds SO 1043 and 
SO 1034 in Tables ^ arid EJ) and allowing for a large 5 % in the rest (almost unconstraining 
this subset). Doing so we observe in (bl) and (b2) curves that are still negative in the region 
where the nominal cash flows of the 0-spread bonds take place. 



9 



0.010 
0.009 
0.008 
0.007 
0.006 
0.005 
0.004 
0.003 
0.002 
0.001 
0.000 



0.006 
0.004 
0.002 
0.000 
-0.002 
-0.004 
-0.006 



spread =0 | (al) 




-■-SO 1033 
-•-SO 1042 
-A- SO 1035 
-T-SO 1044 
-♦- SO 1038 
— «-SO 1037 
SO 1040 
-•-SO 1043 
-+- SO 1034 



5 10 15 20 

Vy/5 (year"') 
spread =0] (bl) 



mm 



=*=*= 



=1=1 



fy75 (year" 1 ) 



0.005 
0.004 
0.003 
0.002 
0.001 
0.000 



I spread =0.0051 (a2) 

mi i 



| spread =0.01] (a3) 



1 - 0.003- i * * — * 







★ — * E 



♦ 
-4 

-! 



0.005 
0.004- 
0.003 
0.002 
0.001 
0.000 
-0.001 
-0.002 
-0.003 
-0.004- 
-0.005 



1 2 3 

i¥f (year" 1 ) 
spread =0.0051 (b2) 
: **** -*-* — * — * — * — 



W (year" 1 ) 



0.002 
0.001- 
0.000 



0.005 
0.004 
0.003 
0.002 
0.001 
0.000 
-0.001 
-0.002 
-0.003 
-0.004 
-0.005 



• . • • 



2 3 

"iW (year 1 ) 



spread =0.01 (b3) 



< 



2 3 

i¥v (year" 1 ) 



Figure 4: Plots indicating some statistical features of price prediction errors as a function 
of y / 7/V- These prediction errors are obtained removing the indicated bonds from the con- 
straining dataset. The statistical sample comprise the forward curves corresponding to 335 
consecutive trading days starting on Wednesday, November 08, 2000 and ending on Thurs- 
day, March 07, 2002. Within this set of days we calculate the relative errors of the predicted 

prices given by (pP red - - pmarket^ /pmarket where pmarket j g the actual price of the bond 

given by the market and pP red - is the price predicted using the prices of all other bonds in 
the dataset. Hence in plots (al), (a2) and (a3) we show the average of the absolute value 
of these errors and the plots (bl), (b2) and (b3) just their average. Consecutive columns 
show averages obtained using spreads of 0%, 0.5% and 1% in the constraining dataset. The 
set of predicted bonds shown in this figure is given by the first 9 bonds of Table [2] . The 
remaining bonds of this table (the last two bonds with the larger maturities) are analyzed in 
Fig.©. Note how prediction errors tend to decrease when spreads are considered albeit not 
significantly (compare with Fig.©). Note also that in the absence of spreads the DF measure 
(v7/V = 00 ) systematically exhibits a better performance than the AD one (yt/V = 0))- 
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Figure 5: In this figure we consider the same kind of statistical features as in Fig.© but now 
the plots correspond to the last two bonds of Table |2] that are the ones with larger maturities 
within our dataset (we have also included the bond with the shortest maturity to facilitate 
the comparison with Fig.©). For the bond with the largest maturity the forward curve has 
to be extrapolated to reach its last cash flows (see Fig.®). We consider two extrapolation 
possibilities, one where we continue the forward rate curve from the last cash flow as a constant 
and the other where the curve is dictated by the minimization of functional W even beyond 
the last cash flow in the constraining dataset. Note that constant extrapolation gives better 
results than W-generated extrapolation in the region close to the AD measure {y/j/tp = 0). 
Again like in Fig.© we observe that in the absence of spreads the DF measure (a/7/V = 00 ) 
systematically exhibits a better performance than the AD one. The most striking difference 
with Fig.© is that here we observe that for these two long maturing bonds the introduction 
of spreads drastically reduce the prediction error all along the family of measures. 
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Since the optimization problem we are solving is non-convex we can not discard the 
presence of several local minima (this is independent of the presence of the positivity 
constraint JE]). However, we have tested our algorithm starting from different seeds 
and in all cases we have arrived at the same minima. These tests included hundreds 
of searches with initial log-prices given by p^ = p b j + x ^p" — p^j (with x € [0, 1] a 
flat random variable) and initial forward rates given by several constant and oscillating 
functions. With this comment we only want to convey our practical experience and by 
no means we intend to say that we have exhaustively explored the presence or absence 
of local minima. 

It is clear that one way to avoid negative rates is just by increasing spreads by hand. 
The advantage of using real spreads at a given moment is that one can be sure of 
being consistent with market prices. If one observes in Fig.(J2J) the large forward rate 
variations taking place for different spread patterns no doubts should remain about the 
relevance of a careful treatment of this issue. 

From Fig.© we conclude that the inclusion of spreads can remarkably improve the 
accuracy of resulting forward curves in the prediction of market prices of long maturing 
bonds. In JE] it was pointed out that the inclusion of spreads notably improved the 
smoothness of the forward curve. To the knowledge of the authors this is the first 
time it is shown that their presence also improves the prediction accuracy. For that 
reason we believe that spreads should be considered even when the market data does 
not provide such information. In that case the approach should consist in using a 
cross validation technique to asses the optimal spread minimizing an error criterion 
based in the prediction of market prices (for cross validation methods see for example 
|2U1 I21|). One possibility is using the well know "leave-one-out" cross-validation to 
select the optimal spread much in the spirit suggested by Figs.® and ©• Given a 
constraining set of k products, the method actually consists in obtaining k forward 
rate curves for a given spread, each curve leaving out one of the constraining products 

(bonds in our case) but using only the k omitted products to compute an error criterion 

k 

like \pf red ' - P^ aTket /p™rket_ Thus we obtain a quantitative criterion to select 

i 

an optimal level of spread 6 . For example for our dataset and our family of models 
(measures) it can be conjectured from Figs.© and © that this optimal level of spread 
is typically around 0.5%. 



This technique is stated here only to illustrate the existence of an optimal non vanishing level of spread 
whenever no real spreads are provided by the market. We do not intend to provide an efficient procedure to 
find this level. 
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Figs.® and © strongly suggest that for low spread patterns DF measure is more 
accurate than AD one. For bigger spreads results do not clearly favor any particular 
measure. 
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Appendix A. Newton steps 

Since the objective function Z in Eq.@ is a non-quadratic function of f r , pj we will use 
an iterative quadratic approximation (Newton steps) to find its minima. We start with 
feasible seeds ff^ and p^ fulfilling inequality constraints (@J and we set up initial log-barriers 
coefficients (J- } > and pP > 0. For newton iteration number s we define 



r ' 



pj -Pj +a j ' ( 8 ) 

with s = 1,2, • • • . The hat over / and p indicate that at each step s, /W and p$ are the 
minima of the quadratic approximation and may not fulfill the inequality constraints Q • To 
assure constraints (j3J) are fulfilled a final redefinition — ► /W, p^ — > p^ is necessary after 
each Newton step. This redefinition is explained in appendix C. Expanding Z up to second 
order in and cr[ s ' we write 



1 ■ ' 0(2) 

= l A is]T Q [s] A [s] + A [s]T B [s] x [s] + A [s]T c [s] + x [s]T a [s] 

+ a MT X [s] + }_ a is]T M ls] a ls] + a ls]T D [s] + (9) 

where collects all terms not depending on A^, o"M or AW. We will use square brackets 
around Newton step indices and parenthesis around dynamic programming ones. Let us now 
work-out the matrices involved in the quadratic approximation. From Eqs.0 and Q we 
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immediately obtain 
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Ai S| 1 / A 



r=l \ Jr 

m 



+ o (aW 3 ) 



yW = _^] E ( ln (^_^) +ln (^_^)) 



j=l \Pj 



Pj ''PI 



~fs — 11 m 

i=i 



J? + o(.f) 
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X(r,x,y) , . • <■,., 



Thus denning 

1 x < r < y f 1 i = j 

, h.j ■= \ 
otherwise [ otherwise 

from above expansions and Eq.@ we immediately obtain 

nj -l ( R ^]- 1 



= irX (r, 1, Rtf] ~ 1) - -r-fi a v ex P 



[s] [s-l] 



r=l 



(1 — <5r,l) (1 - <^r,2) £r-l . (1 ~ <5r,n) (1 ~ <5 r ,i) £ r . (1 ~ <5r,n-l) (1 ~ ^r,n) £r+l 



(tr-2 + Zr-lf 



+ 



f 2 F 2 



+ 



(^+er+l) 2 e 2 



£r+l<V+2,a: &C+1 <V— 2,a 



(1 — ^ r ,i) ^jr+M 

(£r + Cr+l) 2 £r+l£r (6r + 6r+l) 2 £a;+l£z (Cr-1 + £r) £r-l£r 
(1 — ftr.l) Cx^r-l,x _ (1 — ^r,n) £rjr-l£: _ (1 ~ &x,n) Cx^ r +l,x 
' (6r-l + 6r) 6r-l^ ~ (&•-! + £r) £r£ r 2 _i " (6.-1 + £r) 6r^_l 

+ 7, 



(1 — $r,n) 72^r,x + (1 — <5r,l) 75 ~~ Tf^r+l.x — 7^7— l,:r 



>-l] 



+ 7" 



r-l 



£r-l 



f [ S -l] _ f [s-l] n r \ 

4r-l - 7 7 7 & 



' /i s_1] - , - fti 1] \ (i - M (i - ^ 



£r-l 



+ 



fr-2 / (£r-2 + &-l) fr-l 



Cr-l 



, , / /r+2 ^ ~~ /r+1 ^ . fr ~ /r+1 \ 0- ~ ^r,n) (1 — Sr,n-l) <- 

+ 4 ¥> ; 1 7 — ; — ttz Cr+i 



(£r + £r+l) £r 



4 ^ ( fl7i ] -f^ 1] + fti 1] -fl s - 1] \ (i-M(i-M Cr> 



r-l 



1 - r,t«-l] 



1 



1 



+ 



s-l] 



+ 



6 3,k 



[s-l] 



Appendix B. Dynamic Programming for a quadratic objective 



Given positive definite symmetric matrices Q G R n ' n and M € R c,c and the n-vector C and 
c- vector D, we want to minimize the objective function 

W (A, a) := ^A T QA + A T C + X -o T Ma + a T D, 
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subject to the set of constraints 

n 

a i + a i + ^2^kB k ,j = 0, j = l,...,c, 

k=0 

thus defining a new objective function 

Z (A, a, A) := W (A, a) + (A T B + a T + a T ) A. 



(13) 



For any Q we define 



'ifc-l 



d := max 



<S(s) 



1 x + 
x = 



where e.g. (f = 0, 1 for diagonal, tridiagonal matrices respectively. When d « n and c « n 
then Z can be solved efficiently with dynamic programming. To set up the notation and for 
those readers not familiar with dynamic programming let us here briefly explain the basics 
of this well known method JH]- We start defining 



:= Z, 
gig-i) .— gig) 



A 9 =A* 



? = 1, 



where A* satisfies 



Va„2 (?) 



A 9 =AJ 



0. 



,n, 



,n. 



(14) 



(15) 



For q = 0, . . . , n we use the inductive hypothesis 

= ^A T Q^A + A T B^X + A T C^ + X T a^ 
+ ^A T G (,?) A + b {q) + a T A + ia T MCT + a T D, 



(16) 



where G G i? c ' 9 and G i? c ' c . The final step of this backwards process 

consists in obtaining a* and A* satisfying 



From Eqs.p5j) and (|16fl we obtain 

1 



9-1 



A=A* 



(17) 



A* 



n(g) I £ QgAr + ^A + ^l' 



(18) 



16 



and plugging Eq.lfTHjl into Eq. lfTHl we obtain 

1 1 



^■(9-1) _ g;{q) 



A 9 =o 2 q(i) 



1 i- 1 ( n (q) n 

{r,s}=l \ "^9,9 



9-1 

9 ' s 0(7-- g + d)0(s- 9 + d) A a 



8- 1 // /)(?)n(!) \ n <,l) r {l) 

r=l \ \ V9,9 / <A,9 



+ i A T f G<9) - ^ g ^ g I A + f a^ T - 



1 r* v 

" 1 A--^ + 



2 I nil) j \ n (q) j 2 r,^) 

^q,q / \ ^9:9 / V' 



(/,(/ 



+ a T \ + ^a T Ma + cr T D, 



hence obtaining 



nf?) r(«) 

flfe" 1 ) = flW - (r - , + d) , 

Q?,9 



WT 1 n («) 



B 



( 9 -i) = a (g) _ WJZi_ 



Wq,q 



b {q-l) = 5(9) 1 C 9 



(9)2 



where 



^9,9 
9(x) :=i 



1 x > 0, 
x < 0. 



Now from Ea.(|17|) we obtain 



Ma* + A* + D = 0, 
G (o) A * + a * + (o) =0 
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or 



a* = —M~ x (A* + D) , (25) 
and using Eq. (|22|) we obtain 



y g =i y y g =i vg,<? y 

Finally from A* and Eq.(|25|) we obtain a* and then using Eqs. (jl8j) and (jl9!21j) we obtain A* 
moving forward in the q index. 

B.l Computing times of a dynamic programming iteration 

The time necessary to compute all QW j s proportional to 

r, = ( n - 4 ^±i> + v iii±i) = ^±i> In - 2 -i±l). 

2. 2 2 \ 3 / 

i=i v y 

To calculate the inverse of M~ x - that is a c x c symmetric matrix, we require a 

computing time proportional to 

T2 = l(c 3 -c). 

All CW.BW and a® require, respectively, computing times proportional to 

d-l , 

T 3 = (n - d) d + = - (2n - d - 1) , 

i=l 

T 4 = cT 3 , 
T 5 = era. 

Finally the calculation of G^ requires a computing time proportional to 

In the forward rate calculation typically we have n > c > d (in Eq.Q we have d = 2) and 
therefore the maximum delay would be given by Tq. 

Appendix C. Fulfilling inequalities and updating log-barriers 

The Newton step obtained from Eq.(jHJ may not fulfill the inequality constraints (jlj). To 
satisfy such constraints we search for a aM in the interval [0, 1] such that 

ft 1] + « [S] AW > 0, p) < p [ ;- 1] + a^af < p]. (26) 
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In order to do this we first determine the maximum a^Lx m the interval [0, 1] satisfying 
Ea, (|26|) . That is, given the sets of points 



we have 



-4; 



[s-iy 



A := {r, /W < o} , 



r£A 



A 



, mm 



' b [b-1] 
Pi ~ Pi 



, mm 

jeA> 



Pj ~ Pj 



and then we take 



with < (3 < 1 (in our implementation we have taken (3 = .9 that is a standard election in 
the optimization literature). Once we have we define 



W [s-i] , [si W 
//W = max(f (aWj/il'- 1 ^), 

/ u W =max(^(aH) / u [s " 1 Umin), (27) 

where /x m i n and /2 m i n are positive small values guaranteeing that matrices Q and M are 
positive definite and \& : [0, 1] — > i? + is a monotonically decreasing function satisfying ^ (0) = 
1 . In our implementation we have taken ^ of the form 



f[s-l 



^(a) = (l-/)(l-a) f + /, 



(28) 



with Z = 10 2 , £ = 1. The value of / controls how fast barriers are reduced. In our imple- 
mentation we have found that adequate values for I range 10" 3 < I < 10 -1 . The value of £ 
controls the non-linearity of ^f. We have found that the simple linear response provides good 
performance albeit convergence time is not significantly affected for £ in the range 0.6 < £ < 2 
In this way we iterate the algorithm nu times until a given termination criterion is met. 
Defining 



4iw :=l*(w(fW))-ln(w(fW 

( <-i 

pf - ln {vf) + E # ] s> 



max 
3 



r=l 
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we have chosen the following termination criterion 



THt = N max or | 
and 



< W zero or ( 



< 0\ U W ancL °lnVK 



- emax 
<- °\nW . 



nu > N min and fi^ < 



In our 



implementation we have taken 



/i max and /2 [s] < /i max and e W < e 



max} • 



(29) 





= 10" 


> /^min — 10 


10 

5 


Mmax — 


io- 


AM 


= 10+ 


> /^min — 10 


10 

5 


Mmax — 


10" 


^max 


= 10" 


8 srmax _ i «- 


-2 

5 


zero — 


10" 


A'min 


= 5, 


Amax = 60. 









In IF and ^Lw' 



(30) 

Obviously there is considerable latitude to change the heuristic values assigned to the 
above parameters. Let us finish this appendix making some comments regarding their ro- 
bustness. 

iV m i n is there to guarantee a minimum number of iterations so as to have 6 } 
well defined and to avoid premature termination in the improbable case where the other 
criteria incorrectly suggest convergence. For this purpose is enough to take 2 < N m i n < 5. 

Amax serves as a maximum limit to secure termination even if convergence is not achieved 
and therefore an alarm should be provided whenever nu = N max . From our experience we 
observe that is more than enough to take 50 < iV max < 100. 

W zero sets our precision to consider a given forward rate curve as a straight line. The 
order of magnitude of W zero should be taken much lower than the typical order of magnitude 
of the observed optimal W. The value of the optimal W depends not only on the constraining 



data but also on 7 and tp. We have adopted the practise of spanning the range < yT/y? < 1 
year -1 taking < 7 < 1 year 3 , tp = 1 year 5 and the range 1 year -1 < yT/V < 00 taking 
7 = 1 year 3 , < cp < 1 year 5 . With this convention we have found reasonable to take 



10- 10 < W zero < 10" 8 for the full range < ySV < 00. 

^hiW se ^ s ^ e max i mum variation of In (W) between newton steps that is accepted before 
termination. Note that in (|29|) to have ( 5^j w < S^yy and <5 ml y < is a necessary but 

not sufficient condition for termination. Therefore sending — > 00 has no major impact 

and is equivalent to rely only on the error eM and the barrier coefficients as indicators of 
convergence. On the contrary excessively reducing can generate unnecessary iterations. 
We have tested that for values satisfying Sf^y > 10 -5 we do not have any drastic increase in 
convergence time. 

e max controls the error in the constraints and is the most important parameter in (|29fl . 
A too large value of e max reduces the accuracy of the result and a too small value can give 
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rise to unnecessary iterations. We have observed acceptable results for e max in the range 

in -4 < e < irr 10 

lu fc max lu 

/^max and /i max control the maximum allowed values for the logarithmic barriers imple- 
menting the positivity and spread constraints respectively. Large values for these parameters 
can make log barriers to have a residual influence in the feasible region. Acceptable values of 
these maximum weights range in fi m i n < /i m ax ^ 10~ 6 and /2 m i n < /J m ax < 1CP 6 . As explained 
above keeping parameters /i m i n and /i m i n positive guarantees that matrices Q and M are pos- 
itive definite which in turn is sufficient to guarantee the existence of an optimal solution in 
each iteration. Hence /i m i n and /2 m i n should be chosen as small as possible without interfering 
with the numerical stability of the algorithm. Using double precision in our program we have 
found the values given in l|29|) as a good compromise. 

Finally /^°] and jffl are the initial values for the log barriers coefficients. Taking very 
large values for /jM and ft,^ increases the convergence time because we need more time to 
reduce the barriers. Taking too small values for /j^ and /ii ' also increases the convergence 
time because time is wasted exploring unfeasible solutions. Moreover, we have observed that 
convergence and stability is improved if the contributions to the objective of the two barrier 
terms are kept balanced. This is achieved setting 

£[0] „ n_^[0] ( see Eqs.(HIJ and (d) and 
using the same updating factor ^ (aM) for /jM anc [ pj?] ( see Ea.(|27|>). Keeping ~ ^mA 1 ^ 
we have found that the time of convergence is stable for in the range 10 4 < fJ- ^ < 10 +1 . 



Appendix D. Bond tables 

In this work we have used the following data tables and conventions. 
The price of the bond j is calculated using the formula 

i=i l i + TooJ 

where Nj is the nominal amount (SEK 40 millions for all bonds in Table Cj is the coupon 
rate (given in TableEJ, rij is the total number of remaining coupons (each paid at time T-), r,- 
is the quoted rate given in Tableland ATj is a time difference between T- and the settlement 
day. This time difference is calculated according to the ISMA 30E/360 convention defined 
as follows: given two dates (rfi,mi,yi) and (cfo, 7772, 7/2), their ISMA 30E/360 time difference 
AT is given by 

if d\ = 31 set d\ to 30, 

if ^2 = 31 set d.2 to 30, 

30 (771-2 — 777i ) + di — d\ 



AT = y 2 - y\ + 



360 
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Date\Bond (SO) 


1033 


1042 


1035 


1044 


1038 


1037 


1040 


1043 


1034 


1045 


1041 


Friday 06 


4.86 


4.92 


5.06 


5.15 


5.26 


5.27 


5.355 


5.4 


5.395 


5.46 


5.655 


Monday 09 


4.905 


4.965 


5.11 


5.2 


5.05 


5.325 


5.4 


5.455 


5.23 


5.51 


5.7 


Tuesday 10 


4.885 


4.945 


5.095 


5.185 


4.93 


5.295 


5.395 


5.435 


5.24 


5.495 


5.68 


Wednesday 11 


4.865 


4.92 


5.06 


5.15 


4.93 


5.275 


5.355 


5.41 


5.415 


5.465 


5.65 


Thursday 12 


4.835 


4.885 


5.025 


5.125 


4.93 


5.25 


5.355 


5.405 


5.405 


5.465 


5.66 


Friday 13 


4.84 


4.89 


5.045 


5. 145 


4.93 


5.27 


5.38 


5.43 


5.43 


5.495 


5.69 


Monday 16 


4.825 


4.87 


5.035 


5.13 


4.93 


5.25 


5.36 


5.415 


5.41 


5.475 


5.66 


Tu c s d ay 1 7 
























Wednesday 18 


4.8 


4.86 


5.005 


5.1 


4.93 


5.22 


5.335 


5.38 


5.395 


5.445 


5.64 


Thursday 19 


4.77 


4.83 


4.97 


5.06 


4.93 


5.165 


5.27 


5.34 


5.34 


5.39 


5.585 



Table 1: Quoted rates for the eleven Swedish Government Bonds used in the calculation of 
the forward rate curves corresponding to Figs. (|ll3j) . All days are from a two- week period of 
July 2001. The description of each of the bonds is given in Table |2j The prices of the bonds 
are obtained from this rates and the data of Table [21 using Ea. (|31j) . 
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Maturity 
(dd/mm/ yyyy) 


annual 
coupon (%) 


□ U 1033 


nr /fir /onno 
Uo/0o/ zU03 


10. 25 


SO 1042 


15/01/2004 


5 


SO 1035 


09/02/2005 


6 


SO 1044 


20/04/2006 


3.5 
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25/10/2006 


6.5 


bU 103f 


Id/Oo/ zUO i 


8 
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nr /Ar /nriAO 
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bU 104 J 


OO /A1 /flAAfl 


5 
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05/05/2014 
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